home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / wilcoxon.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  214 lines
  1. ; $Id: wilcoxon.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7. pro Wilcoxon,Data, Rank,Prob,Names=Names,List_Name=Ln,    $
  8.                   Missing=M,Tail =T,NoPrint=NP 
  9. ;+ 
  10. ; NAME:
  11. ;    WILCOXON
  12. ;
  13. ; PURPOSE:
  14. ;    Test the null hypothesis that two populations have the same 
  15. ;    distribution -i.e. F(x) = G(x) against the alternative that their 
  16. ;    distributions differ in location- i.e F(x) = G(x+a).
  17. ;    Wilcox pairwise tests the populations in Data.  It is preferable 
  18. ;    to the sign_test for non-binary data, since it uses the ranking of
  19. ;    data differences and not just the sign.
  20. ;
  21. ; CATEGORY:
  22. ;    Statistics.
  23. ;
  24. ; CALLING SEQUENCE:
  25. ;    WILCOXON, Data
  26. ;
  27. ; INPUTS:
  28. ;    Data:    Two-dimensional array. Data(i,j) = the jth observation from 
  29. ;        the ith population.
  30. ;
  31. ; KEYWORDS:  
  32. ;    NAMES:    Vector of user supplied names for the populations to be 
  33. ;        used in the output.
  34. ;
  35. ;    LIST_NAME:    Name of output file.  The default is to the screen.
  36. ;
  37. ;      MISSING:    Value used as a place holder for missing data.
  38. ;        Pairwise handling of missing data.
  39. ;
  40. ;    TAIL:    If set, tail specifies the type of test.  If Tail = 1 then a 
  41. ;        one-tailed test is performed.  Likewise, Tail = 2 specifies 
  42. ;        a two-tailed test.  The default is the two-tailed test.
  43. ;
  44. ;      NOPRINT:    A flag, if set, to suppress output to the screen or a file.
  45. ;           
  46. ; OUTPUT:
  47. ;    Table written to the screen showing for each pair of populations the
  48. ;    total of the rankings.  Also, table of probabilities for each 
  49. ;    population pair giving the significance of the results in the first
  50. ;    table. These probabilities are based on normal approximations.  For
  51. ;    small sample sizes (less than 26), consult a table for the Wilcox 
  52. ;    Test in any general text on statistics.
  53. ;
  54. ; OPTIONAL OUTPUT PARAMETERS:
  55. ;    Rank:    Two-dimensional array of ranking sums.
  56. ;        Rank(i,j) = sum of ranks of population i where observations
  57. ;        in i are ranked with those in population j.
  58. ;
  59. ;    Prob:    Two-dimensional array. Prob(i,j) = probability of Rank(i,j) 
  60. ;        or something more extreme. -1 is returned for small 
  61. ;        (less than 26) populations.
  62. ;                          
  63. ; RESTRICTIONS: 
  64. ;    All populations have the same sample size.
  65. ;
  66. ; COMMON BLOCKS:
  67. ;    None.
  68. ;
  69. ; PROCEDURE:
  70. ;    For each pair of populations,  differences between corresponding
  71. ;    observations are computed and the absolute values of the differences
  72. ;    are ranked.  Ranks of positive differences are summed and likewise,
  73. ;    those of the negative ones.  The probability of the sums is computed
  74. ;    under the assumption that the distributions are the same.  These 
  75. ;    probabilities are approximated with an approximately normal test
  76. ;    statistic for populations greater than 25.
  77. ;-
  78.  
  79.  
  80.  
  81.  
  82. On_Error,2
  83. SD= size(Data)
  84.  
  85. if( N_Elements( Ln) NE 0) THEN openw,unit,/get,Ln else unit=-1
  86.  
  87. if ( SD(0) NE 2) THEN BEGIN
  88.    printf,unit, 'wilcoxon- Data array has wrong dimension'
  89.    goto, DONE
  90. ENDIF
  91.  
  92.  if(N_Elements(T) EQ 0) THEN T=2   $
  93.  else if (T GT 2) OR  (T LT 1) $
  94.       THEN T =2
  95.  
  96. C=SD(1)
  97. R= SD(2)
  98.  
  99. Rank = Fltarr(C,C) 
  100. Prob =  Rank +1.0
  101. for  i =1l,C/2 DO  BEGIN             ;cycle through populations
  102.  
  103.    D1 = Shift(Data,i,0)
  104.    Temp = Data - D1                  ; compute differences
  105.  
  106.    if (N_Elements(M) NE 0) THEN BEGIN  ; Handle missing data
  107.      here = where(Data EQ M,count)
  108.     if count  NE 0 THEN    BEGIN
  109.       Temp( here) = 0
  110.       Temp( where(D1 EQ M)) = 0
  111.     ENDIF
  112.    ENDIF
  113.    
  114.  
  115.    for j = 0l, C-1 DO BEGIN             ;compute rankings
  116.      k = (j + (C-i)) mod C
  117.      accept =0
  118.      SortT =Temp(j, sort(abs(Temp(j,*))))
  119.                                     ;sort absolute values
  120.                                          ; of differences
  121.      pos = where(SortT NE 0,countp)    ; discard equal
  122.                                        ; observation
  123.  
  124.     if countp eq 0 THEN BEGIN
  125.          Rank(k,j) = 0 & Rank(j,k) = 0
  126.          prob(k,j)=  1 & prob(j,k) =1
  127.     ENDIF else BEGIN
  128.         SortT = SortT(pos)
  129.                                 ; if all equal, accept NUll 
  130.                     ; hypotheses
  131.         PopSize = N_Elements(SortT)    
  132.         Pos =float( where(sortT GT 0,countp))
  133.                                    ; get rankings of positive
  134.         Neg =float( where (sortT LT 0,countn))
  135.                                    ; of negative
  136.         Repeats =where ( abs(sortT) EQ abs(shift(sortT,-1)),count)
  137.                                    ;look for repetitions to average
  138.  
  139.         if count NE 0 THEN BEGIN
  140.           here = 0
  141.  
  142.           while ( here LT count-1) DO BEGIN
  143.             diff = abs(SortT(Repeats(here)))
  144.             set1= where( abs(SortT) EQ diff,count2)
  145.             if countp NE 0 THEN    $
  146.               posr = where(SortT(pos) EQ diff,count3)  $
  147.                          ELSE count3 = 0
  148.             ave = repeats(here) + (count2-1.0)/2.0
  149.             if count3 NE 0 THEN Pos(posr) = ave
  150.  
  151.             if countn NE 0 THEN posr =            $
  152.                     where (SortT(neg) EQ -diff,count3) $
  153.               ELSE count3 = 0
  154.  
  155.             if count3 NE 0 THEN Neg(posr) = ave
  156.             here = here + count2 -1
  157.  
  158.        ENDWHILE
  159.       ENDIF
  160.  
  161.       if countp NE 0 THEN  Rank(j,k) =      $
  162.              Total(Pos+1) ELSE RANK(j,k) = 0
  163.       if countn NE 0 THEN  Rank(k,j) =      $
  164.                       Total(Neg+1) ELSE RANK(k,j) = 0
  165.        Z = (Rank(k,j) - (PopSize * (PopSize+1)/4.))/       $
  166.            sqrt(popsize*(popsize+1.)*(2*popsize+1.)/24.)
  167.       if popsize gt 25 THEN  Prob(j,k) =( 1 - Gaussint(abs(Z))) $
  168.          else Prob(j,k) = -1
  169.        if T EQ 2 and popsize gt 25 THEN $
  170.                Prob(j,k) = 2*Prob(j,k)
  171.        Prob(k,j) = Prob(j,k) 
  172.     ENDELSE
  173.    ENDFOR
  174. ENDFOR
  175.  
  176.  SN =Size(Names)
  177.  if (SN(1) EQ 0) THEN BEGIN
  178.      I = INDGEN(C)
  179.      Names=['pop'+StrTrim(I,2)]  
  180.  ENDIF ELSE                  $
  181.    if ( SN(1) LT C) THEN BEGIN
  182.      I = Indgen(C)
  183.      printf,unit,'sign_test- missing names'
  184.      Names=[Names, 'pop'+StrTrim(I(SN(1):C-1),2)]
  185.   ENDIF
  186.  
  187.  If NOT KEYWORD_SET(NP) THEN    BEGIN
  188.   printf,unit, "Wilcoxon"  
  189.   printf,unit, " Table of Rank Sums"
  190.   printf,unit, " "
  191.   printf,unit, format ='(8X,16(A15,2x))', Names
  192.  
  193.   for i= 0l,C-1 do                       $
  194.      printf,unit, format='(A8,16(G15.5,2X))',   $
  195.                           Names(i),Rank(*,i)
  196.  printf,unit, " "
  197.  printf,unit, "Table of Probabilities:"
  198.  if (T EQ 2) THEN       $
  199.    printf,unit, "Two- Tailed" $
  200.  else printf,unit, "One-Tailed"
  201.  
  202.  printf,unit," "
  203.  printf,unit, format ='(8X,16(A15,2x))', Names
  204.  for i= 0,C-1 do                      $
  205.      printf,unit, format='(A8,16(G15.5,2X))',     $
  206.             Names(i),Prob(*,i)
  207.  ENDIF
  208.  
  209.  
  210. DONE:
  211.    if ( unit NE -1) THEN Free_Lun,unit
  212.     RETURN
  213.     END
  214.